!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_variables (ro, mx, my, mz, bx, by, bz, en, gx, gy, gz, po, x, y, z, nx, ny, nz)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz

   real(8), intent(in) :: x(nx), y(ny), z(nz)
   real(8), dimension(nx, ny, nz), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, gz, po

   real(8), dimension(nx, ny, nz) :: pr, te
   real(8), dimension(ny_total_finest) :: te_total_y

   integer(4) :: i_current, i_charge, j_current, j_charge, j_total
   integer(4) :: nx_current, ny_current, nx_charge, ny_charge, zeroy

   real(8) :: amp_b, d_current, d_charge, start_x, start_y, start_z, ri, ri3, r0_x, r0_y, r0_z
   real(8) :: g0, te_chr, te_cor, h_cor, w_tr, dy

!----------------------------------------------------------------------------------------------------------------------------------|
!  1d te, ro distribution
!----------------------------------------------------------------------------------------------------------------------------------|
   g0 = -1.0d0
   te_chr = 0.6d0
   te_cor = 100.0d0
   h_cor = 8.5d0
   w_tr = 1.0d0

   if (.not. allocated (y_total)) allocate (y_total(ny_total_finest))
   if (.not. allocated (ro_total_y)) allocate (ro_total_y(ny_total_finest))

   if (ro_total_y_flag) then

      zeroy = (nint((ny_total - 2 * ng) * zero_posy / rangey) + ng) * 2 ** (n_levels - 1) + 1
      dy = rangey / dble(ny_total - 2 * ng) / 2.0d0 ** (n_levels - 1)

      y_total(zeroy) = dy * 0.5d0

      do j = zeroy + 1, ny_total_finest
         y_total(j) = y_total(j - 1) + dy
      enddo
      do j = zeroy - 1, 1, -1
         y_total(j) = y_total(j + 1) - dy
      enddo

      do j = 1, ny_total_finest
         te_total_y(j) = te_chr + (te_cor - te_chr) / 2.0d0 * (tanh((y_total(j) - h_cor) / w_tr) + 1.0d0 + eps)
      enddo

      ro_total_y(zeroy) = (2.0d0 * (te_total_y(zeroy - 1) + te_total_y(zeroy)) / 2.0d0 + dy * 0.5d0 * g0) /                        &
         (2.0d0 * te_total_y(zeroy) - dy * 0.5d0 * g0) * 1.0d3
      do j = zeroy + 1, ny_total_finest
         ro_total_y(j) = (2.0d0 * te_total_y(j - 1) + (y_total(j) - y_total(j - 1)) * g0) /                                        &
            (2.0d0 * te_total_y(j) - (y_total(j) - y_total(j - 1)) * g0) * ro_total_y(j - 1)
      enddo

      do j = zeroy - 1, 1, -1
         ro_total_y(j) = (2.0d0 * te_total_y(j + 1) - (y_total(j + 1) - y_total(j)) * g0) /                                        &
            (2.0d0 * te_total_y(j) + (y_total(j + 1) - y_total(j)) * g0) * ro_total_y(j + 1)
      enddo

      ro_total_y_flag = .false.

   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  initial value
!----------------------------------------------------------------------------------------------------------------------------------|
   j_total = 1
   do j = 1, ny
   do k = 1, nz
   do i = 1, nx

      gx(i, j, k) = 0.0d0
      gy(i, j, k) = g0
      gz(i, j, k) = 0.0d0

      te(i, j, k) = te_chr + (te_cor - te_chr) / 2.0d0 * (tanh ((y(j) - h_cor) / w_tr) + 1.0d0 + eps)

      do while (y_total(j_total) .lt. y(j))
         j_total = j_total + 1
      enddo

      if (j_total .eq. 1) then
         ro(i, j, k) = ro_total_y(j_total)
      else
         ro(i, j, k) = (ro_total_y(j_total - 1) * (y_total(j_total) - y(j)) +                                                      &
                        ro_total_y(j_total) * (y(j) - y_total(j_total - 1))) /                                                     &
                       (y_total(j_total) - y_total(j_total - 1))
      endif
      pr(i, j, k) = ro(i, j, k) * te(i, j, k)
      mx(i, j, k) = 0.0d0
      my(i, j, k) = 0.0d0
      mz(i, j, k) = 0.0d0
      bx(i, j, k) = 0.0d0
      by(i, j, k) = 0.0d-1
      bz(i, j, k) = 0.0d0

   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  magnetic chage
!----------------------------------------------------------------------------------------------------------------------------------|
   nx_current = 24
   ny_current = 24
   d_current = 100.0d0

   nx_charge = 8
   ny_charge = 8
   d_charge = 6.0d0

   start_x = -(nx_current - 1) * d_current / 2 - (nx_charge - 1) * d_charge / 2
   start_y = -ny_current * d_current / 2 - (ny_charge - 1) * d_charge / 2
   start_z = -3.6d0

   amp_b = 1.8d5 / nx_current / ny_current / nx_charge / ny_charge

!----------------------------------------------------------------------------------------------------------------------------------|
!  canopy
!----------------------------------------------------------------------------------------------------------------------------------|
!   do j_current = 1, ny_current
!   do i_current = 1, nx_current
!   do j_charge = 1, ny_charge
!   do i_charge = 1, nx_charge

!      r0_x = start_x + (i_current - 1) * d_current + (i_charge - 1) * d_charge
!      r0_y = start_y + (j_current - 1) * d_current + (j_charge - 1) * d_charge
!      r0_z = start_z

!      do k = 1, nz
!      do j = 1, ny
!      do i = 1, nx
!         ri = sqrt((x(i) - r0_x) ** 2 + (y(j) - r0_y) ** 2 + (z(k) - r0_z) ** 2)
!         bx(i, j, k) = bx(i, j, k) + amp_b / nx_current / ny_current / nx_charge / ny_charge * (x(i) - r0_x) / ri ** 3
!         by(i, j, k) = by(i, j, k) + amp_b / nx_current / ny_current / nx_charge / ny_charge * (y(j) - r0_y) / ri ** 3
!         bz(i, j, k) = bz(i, j, k) + amp_b / nx_current / ny_current / nx_charge / ny_charge * (z(k) - r0_z) / ri ** 3
!      enddo
!      enddo
!      enddo

!   enddo
!   enddo
!   enddo
!   enddo

   do j_current = 1, ny_current
   do i_current = 1, nx_current
   do j_charge = 1, ny_charge
   do i_charge = 1, nx_charge

      r0_x = start_x + (i_current - 1) * d_current + (i_charge - 1) * d_charge
      r0_y = start_z
      r0_z = start_y + (j_current - 1) * d_current + (j_charge - 1) * d_charge

      do k = 1, nz
      do j = 1, ny
      do i = 1, nx
         ri = sqrt((x(i) - r0_x) ** 2 + (y(j) - r0_y) ** 2 + (z(k) - r0_z) ** 2)
         ri3 = amp_b / ri ** 3
         bx(i, j, k) = bx(i, j, k) + (x(i) - r0_x) * ri3
         by(i, j, k) = by(i, j, k) + (y(j) - r0_y) * ri3
         bz(i, j, k) = bz(i, j, k) + (z(k) - r0_z) * ri3
      enddo
      enddo
      enddo

   enddo
   enddo
   enddo
   enddo

!   do k = 1, nz
!   do j = 1, ny
!   do i = 1, nx
!      ri = sqrt((x(i) - r0_x) ** 2 + (y(j) - r0_y) ** 2 + (z(k) - r0_z) ** 2)
!      ri3 = amp_b / nx_current / ny_current / nx_charge / ny_charge / ri ** 3
!      bx(i, j, k) = bx(i, j, k) + (x(i) - r0_x) * ri3
!      by(i, j, k) = by(i, j, k) + (y(j) - r0_y) * ri3
!      bz(i, j, k) = bz(i, j, k) + (z(k) - r0_z) * ri3
!   enddo
!   enddo
!   enddo

   do k = 1, nz
   do j = 1, ny
   do i = 1, nx

      en(i, j, k) = (mx(i, j, k) * mx(i, j, k) + my(i, j, k) * my(i, j, k) + mz(i, j, k) * mz(i, j, k)) / ro(i, j, k) / 2.0d0 +    &
         (bx(i, j, k) * bx(i, j, k) + by(i, j, k) * by(i, j, k) + bz(i, j, k) * bz(i, j, k)) / 2.0d0 + pr(i, j, k) / gmm1

      po(i, j, k) = 0.0d0

   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
end subroutine model_variables
